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Abstract 

A class of implicit upwind-differencing 
methods for the compressible Navier-Stokes equa- 
tions is described and applied. The methods are 
based on the use of local eigenvalues or wave 
speeds to control spatial differencing of inviscid 
terms and are aimed at increasing the level of 
accuracy and stability achievable in computation. 
Techniques for accelerating the rate of convergence 
to a steady-state solution are also used. Applica- 
tions to inviscid and viscous transonic flows are 
discussed and compared with other methods and exper- 
imental measurements. It is shown that accurate 
and efficient transonic airfoil calculations can be 
made on the Cray-1 computer in less than 2 min. 


Introduction 

The purpose of this paper is to describe a 
class of implicit upwind-differencing methods (IUM) 
for the compressible Euler and Navier-Stokes equa- 
tions. The methods use local eigenvalues or wave 
speeds to control inviscid spatial differencing and 
are closely related to many other recent tech- 
niques. The objective in all these methods is 

to achieve more stable and accurate solutions than 
can be obtained using conventional techniques such 
as the central-differencing method 10 and 
MacCormack's method. 11 In the latter two methods, 
dissipative terms must usually be added to control 
parasitic oscillations, and the choice of the form 
of the dissipation terms and the size of the dissi- 
pation constants is difficult. The more recent 
methods, utilizing local eigenvalues, are naturally 
dissipative and, in principal, do not require the 
addition of extra terms to stabilize calculations. 

In practice, however, these methods also have dif- 
ficulties, especially in regions of the flow where 
the eigenvalues change sign, and again special pro- 
cedures must frequently be introduced to improve 
accuracy. 1-3 

The class of methods reported in this paper 
is an extension of the second-order method reported 
in Ref. 1 where dissipative terms scaled on the 
squares of the local eigenvalues were used. 

Although the method of Ref. 1 is simpler chan the 
present methods, it was found in subsequent appli- 
cations to suffer the same difficulty as that of 
earlier methods, that is, the need for fine tuning 
of dissipation constants to achieve accurate solu- 
tions. Of the class of mechods to be reported in 
this paper, one second-order upwind method (UW2II) 
has been found, through extensive numerical experi- 
mentation, to produce accurate solutions of 
transonic-flow problems without the need for spe- 
cial treatment at normal shock waves or sonic lines. 
For oblique shock waves, the method produces 3mall 
undershoots or overshoots around the shock which 
may be substantially reduced by a simple technique. 
The other mechods reported produce accurate 


solutions, except at normal shock waves, where they 
may be improved by special treatments. 

The methods presented are closely related to 
the flux-vector splitting method of Steger and 
Warming, 3 who used the finite-difference technique. 
The present methods are designed for incorporation 
into the strongly conservative finite-volume tech- 
nique 1 which, in contrast to the former technique, 
retains the property of free-stream maintenance in 
curvilinear coordinates without the need for dif- 
ferencing metrics in the same way as fluxes. 

In this paper, a one-dimensional description 
of the mechods will be made followed by a discus- 
sion of two-dimensional results which include 
inviscid and viscous transonic- flow calculations. 
These are compared with other numerical methods and 
experimental results. A basic conclusion from 
these results is chat accurate and efficient tran- 
sonic airfoil calculations can be made, using the 
full Navier-Stokes equations, in less chan 2 min 
on the Cray-1 computer. 

Theory 

Although the methods will be applied to two- 
dimensional viscous-flow problems in curvilinear 
coordinates, they will be described here in Che 
inviscid, one-dimensional context in order to 
simplify the development. The one-dimensional 
Euler equations may be written 

3 t U + 3 x F - 0 (1) 

where 3 C - 3/3t, 3^ - 3/3x, and 

U - (p,pu,pE) T , F - (pu,pu 2 + p » (pE + p)u) T (2) 

are the conservative state and flux vectors with 
0 * density, u - velocity, p - pc 2 /y * pressure, 
oE - p/(y - 1) + pu 2 /2 ■ total energy, and 
c - sound speed. The nonconservative state vector 
V and the Jacobian of the flux vector A are 
written 

V - (p,u,p) T , A - 3F/3U - R _1 AR O) 

where A - diag(u,u + c, u - c) is the diagonal 
matrix of eigenvalues of A, and R is a similarity 
matrix diagonalizing A. The matrix R may be 
written 

R - QP, P - 3V/3U (*) 



m 

- c 2 

0 
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1 

0 

0 “ 
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0 

pc 

1 
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-u/p 

1/p 
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0 

-pc 

1 
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X 2 

L 2 u 

-Xu 

X 


where X ■ y - 1 and y is the ratio of specific 
heats . 


^Research Scientist. Member AIAA. 

This paper is declared a work of the U.S. Government and 
therefore is in the public domain. 


1 


The sign of A and Che absolute value of A are 
defined by 

S * sgn A - R“ l sgn A R 

(5) 

| A | " abs A - R“ l abs A R ■ SA - AS 

where sgn A and abs A are diagonal matrices con- 
sisting of che signs and Che absolute values of the 
eigenvalues, respectively. 

The basic implicit algorithm, using first- 
order time-differencing, may be written in the 
delta form (leaving spatial differencing arbitrary) 

as 

(I + At3 x A)AU - -At3 x F (6) 

where AU ■ U(x,t + At) - U(x,t) is the delta var- 
iable and At is the time-step. 

In order to completely define che algorithm, 
the spatial derivatives in Eq. (6) must be replaced 
by spatial difference operators. We consider first 
the spatial flux differencing, or the right-hand- 
side term of Eq. (6). It is expressed in the fol- 
lowing generic (finite- volume) form 

3 x F * ^i+i/2 “ ^i-l/2^ /Ax ^ 

where Ax is che mesh spacing and F t+i /2 is ctie 

flux vector defined at che midpoint (or cell face) 
between the mesh points i and i+1. The state 
vector, U(x,t) -*■ U(iAx,nAt) ■ Uj, is defined at 
che mesh points, and the flux vector, is 

defined in terms of U^, U^ +l , U i+2 as 

follows: 

*1+1/2 " 1/2(F i + F i+1 - (8) 

where F i - F(U^) and D^+i /2 la a dissipation 
function to be specified. If D^-h/ 2 " 
scheme reduces to pure central-differencing. We 
note chat the differencing of Eq. (7) is fully con- 
servative, regardless of the definition of D t+i/2* 
since it represents a telescoping sum of terms. 

The dissipation function, Di+ 1/2 D, in 
Eq. (8) will be expressed in terms of spatial- 
difference operators, which are defined as follows. 
We let W be either F, U, or V, and we define the 
operators, 6^, by 

«i w i+i/2 - t * +1/2 (Wi + i - »i> * 


The present methods may be developed in a 
sequence of steps starting with the flux-vector 
splitting method of Steger and Warming. 3 The 
latter method can be expressed in terms of Eqs. (7) 
and (8) , using the dissipation function 

D -^AU +^t 1 |a|u (11) 

and assuming that che flux vectors (F^ and F^+^) in 
Eq. (8) possess the homogeneous property, chat is, 

F - AU. However, it is not necessary to represent 
che flux vectors in this manner and they may be 
used directly. Furthermore, AU and | A | U in 
Eq. (11) may be replaced, using Eq. (5), by F and 
SF so that Eq. (11) can be expressed alternatively 
by 

D-U^F+^SF (12) 

This results in a flux-vector splitting method that 

does not depend on the homogeneous property. 

Unfortunately, che occurrence of the matrix S 
inside che operator JT\ causes difficulties in the 
flux-vector splitting method when che eigenvalues 
change sign. J In this case the differencing becomes 
inconsistent, and special procedures must be used 
to stabilize the calculations. 2 

An alternative form of the method, here called 
Method I, was proposed by Hwang 6 (in one-dimensional 
first-order form for the pseudounsteady Euler 
equations). This form is obtained from Eq . (12) 
by taking che matrix S outside the operator, chat 
is, 

(Method I) D -w*iF + S^F (13) 

where a suitable average must be defined for S. 

This method has the advantage that the differencing 
is more consistent when che eigenvalues change 
sign; numerical experiments in one-dimensic a by 
Huang and in two-dimensions by the present ithor 
using the first-order form and simple (linear) 
averaging produced good, normal shock captures. 
However, the second- and third-order forms, which 
are needed for accuracy in multidimensions, were 
found (by the present author) to be weakly unstable 
at normal shock waves, and special procedures were 
required to remove oscillations. The higher-order 
methods were also tried, using nonlinear averaging, 
and with the same kind of result. This method, 
although che most elegant of the three methods to 
be discussed, was therefore put aside in favor of 
Che second method, now to be discussed. 


^ w i+l/2 * a ^l w i+3/2 “ 5 fc W i-l/2^ 

‘ /f i W i+l/2 * 8<5 i W i+l/2 " Y ^ <5 i W i+3/2 * <5 l W i-i/2^ 


Method II may be obtained from Eqs. (9) and 
(13) by che replacement 

4 1 F * A 5l U - A i+l/2 4lU t+1/2 («> 


* W 


The operators *Jf\F and £^\F are then replaced as 
follows : 


The constants a, 8, and y define Che spatial 
accuracy of the differencing and will be discussed 
later. The subscript Z cakes the values 1, 2, 
and 3, and the matrices T ^ +l /2 are 

by 

T l » I, T 2 * R, T 3 - Q (10) 

A suitable average must be given for T 2 in 
terms of U^, U^ +1 for l - 2, 3, and this averaging 
will be discussed below. 


uTiF - Av4T 1 U - R -1 AR-# : U - R* 1 AU* 2 U 
S^F - SA^U - R~ l |A|R^U - R“ l |A|^ 2 U 


(15) 


so that che dissipation function of Method II can 
be written 


(Method IX) D - R _1 (/Ut> 2 U + |A|./r 2 U) 
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or 


°i+l/2 “ R i+l/2 (A i+l/2^2 U i+l/2 

+ l A l i+ l/2^ U i + l/2 ) 

Method II has been found to produce clean and 
accurate (normal) shock captures in numerical exper- 
iments and is the basis of the results to be 
reported. 

Method III is closely related to Method II but 
uses che nonconservative state vector, V. It is 
obtained through the replacements 

6,U - P~ l 6,V , UfjU - P -1 ur 3 V , ♦ P" 1 ^^ 

(17) 

in Eq. (15). By noting Chat RP -1 • Q, the 
replacements in Sq. (15) can alternatively be taken 
as 

uTjF - A i LJ - R“ 1 AQ - R “ l AUT 3 V 

(18) 

S./KF - sa^u - r- 1 |a|q^v - r- 1 |a|^t 3 v 


so chat che dissipation function of Method III can 
be written 

(Method III) D - R _1 (AuT 3 V + |a|u* 3 V) (19) 

Only limited numerical experimentation has been 
done with this method, and no results will be shown 
for it. 


It is important to note chat che replacements 
in Eqs. (15) and (18), indicated by che arrows, are 
equalities only if A is constant; as a result, 
the methods may be expected to produce different 
resulcs, especially at singular points where che 
eigenvalues change sign. As a further note, we may 
add chat the form of che dissipation functions of 
Methods II and III was inspired by a study of 
Har ten's explicit method 5 ’ 7 so chat a more complete 
understanding of che present method, and perhaps 
even further improvements (e.g., entropy conditions) 
may be achieved by a careful study of his and 
related works. 


In order co completely specify che algorithms, 
che method of averaging, that is, computing 
R i+l/2» e cc., in terms of and Ui+i or 
and Vi+i, muse be given. Although chere are many 
possibilities, we have tried only cwo methods; 
linear averaging using the nonconservacive state 
vector (e.g., u i+l / 2 - (u t + u 1+l )/2 and 
c i+i/2 * ( c i * c i+l)/2) and nonlinear averaging 
after Roe 9 as implemented by Yee et al. 7 In the 
latter case, che replacement 5^F Ad^U becomes 
an exact equality and expresses Roe’s property U. 
In che results co be discussed. Roe's averaging 
has been used. The specific formulas used in com- 
puting and ^i+1/2 f° r Method II are 

given below. 


a 


(p i+1 /0.) 


1/2 


U i+l/2 

°i+l/2 


(au i+1 + u.)/(l + a) 

(y - I) t(aH i+i + H 1 )/(l +• a) 


2 , 
J 1+1/2 J 


( 20 ) 


H - E + p/o 


These equations are used in computing R-i+i/ 2 from 
Eq. (4) where it may be noted chat in che product 
R ■ QP, che density cancels so that an average for 
0 is not required. 

As stated earlier, the coefficients, a, 8, and 
y, appearing in the dissipation operators and 
J* define the spatial accuracy of the flux differ- 
encing. Several options for these coefficients 
are listed in Table 1. 


Table 1 Scheme options and nodal clusters 


Scheme 

Desig- 

Coefficients 

Nodal 

cluster. 

? i+l/2 

nator 

a 

8 

Y 

i - 1 

i i+ 1 

i+ 2 

First- 

order 

upwind 

UWl 

0 

1 

0 

• 

o 

• 

Second- 

order 

upwind 

UW2 

1/2 

1 

1/2 

0 

0 

• 

Third- 

order 

upwind 

UW3 

1/6 

1/3 

1/6 

0 

0 o 

• 

Second- 

order 

dissi- 

pative 

CD2 

0 

e 

e/2 

O 

o 0 

0 


The nodal clusters are indicated for the 
midpoint flux vector, F i+1 / 2 where it is assumed 
chat all che eigenvalues of A^ +l / 2 * are positive. 
For the second-order dissipative scheme, chat is, 
CD2, e is a free parameter.. Assuming 
A - constant, spatial accuracy may be checked by 
Taylor series expansion of Eq. (7) utilizing 
Eqs. (8), (9), and (13). The above schemes may 
be used with any of the three basic methods, which 
will be identified by placing a Roman numeral after 
che scheme designator. For example, the second 
order upwind method II will be denoted by UW2II. 


As mentioned in the Introduction, the capture 
of oblique shock waves is degraded to some extent 
by che appearance of oscillations. These may be 
reduced substantially by a high resolution tech- 
nique in which a switch co first-order, chat is, 
forcing, a, y 0 and 8 •+■ 1, at points of local 
maxima or minima in che pressure. For che second- 
order upwind method this is achieved through the 
following formulas: 


6'P 

6"p 


6p 


i+l/2 


Y - max [0.5, 0.5 - (6'p - 6 ,f p) / g ' ] 

+ l <Sp i+l/2 1 + |6f W 2 l 

l <Sp i-l/2 + aP i+l/2 + 5p i+3/ 2 l 


The parameter e' is a small constant 
(e 1 (0.01 - 0.001) p,,,) which controls che rate of 

approach co the first-order relations at those 
points where che pressure departs from monoconicity , 
that is, 6’p > <5"p. The designator HR will be 
appended to che basic designator when chis option 
is used, for example, (JWIIHR. 
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The spatial differencing of the implicit or 
left-hand-side terms of Eq. (6) follows the devel- 
opment given in Ref. 1. First-order upwind differ- 
encing in an (approximate) nonconservative scalar 
tridiagonal form is used. The algorithm may be 
written 


R _1 [I + At(A + V v + A - A X )]RAU - -At3 F (21) 

X A X 

where 


A” ■ (At | A | ) / 2 

V W - (W i - W i-L )/Ax 

*X W - (w i + l - W i )/4x 

Advantages of this formulation are 1) it sub- 
stantially reduces computing time compared with the 
more exact block-tridiagonal form, 2) the upwind 
differencing is dissipative which, coupled with 
the upwind dissipative differencing of the right- 
hand side, enhances the (linear) stability of the 
overall algorithm, and 3) viscous terms may be 
included by means of a simple technique. The main 
limitation of the formulation is that time accuracy 
is lost to some extent, but this is of no conse- 
quence in steady-state problems, which are the 
principal aim of this paper. 

The inclusion of viscous terms is described 
as follows. The one-dimensional Navier-Stokes 
equations may be written 


where 7 X and A x are defined by Eq. (21). Conven- 
tional second-order central differencing is used 
for the viscous-flux differencing on the right-hand 
side of Eq. (24). 


The numerical method, as outlined in the pre- 
ceding paragraphs, may be generalized to multi- 
dimensional curvilinear coordinates utilizing the 
finite-volume technique, following the development 
of Ref. 1. Boundary conditions based on the method 
of characteristics are also described in Ref. 1. 


Spatially varying time-steps are used to 
accelerate iterative convergence to a steady-state 
solution. The procedure used here is a modifica- 
tion of that described in Ref. 1. Expressed in 
terms of two-dimensional Cartesian coordinates, 
che relationship used for the local time-step is 


At - CFL 


min 


(■ 


Ax 

|u| + c 



* CFL 


At 


local 


where CFL is che local Courant number, usually 
taken to be 0(5-10) for most inviscid problems. 

For viscous- flow problems at high Reynolds numbers, 
where a finely spaced mesh is used adjacent to 
solid surfaces, the above formula leads to very 
small values of At in the neighborhood of the 
surface. In these cases, the formula is modified 
by placing a lower limit on At, chat is. 

At - CFL • max(At locai , Ato) 


3 U + 3 (F - B3 (/) 
t x x 


( 22 ) 




-u u 


0 

u* 


Ju" - u')u 2 - u”e (u' - u")u u'l 


(4/3)u , 


(y/Pr)u , 


Pr • C u/K 
P 


where u is che molecular viscosity, Pr is che 
Prandtl number, and e is the specific internal 
energy, e - c 2 /y(y - 1). The viscous analog of 
che time-differenced delta-form algorithm, that is, 
Eq. (6), becomes 

[I ^ At3 x (A - B3 x )]AU - -At3 x (F - B3 x U) (23) 


where on che left-hand side of this equation the 
matrix B, like A, is evaluated at che time-level n. 
The approximate scalar-tridiagonal form of this 
equation is derived by replacing B by che matrix 
vl, where I is che identity matrix, v - Umav/o , 
and u max the largest eigenvalue of B, chat is, 

"‘max * nax(w , l p") . In this case, che above equa- 
tion, expressed In nonconservacive form, becomes 

R _1 [I + At(A3 x - v3*)]RAU - -Aca^F - B3 X U) (24) 


where Atg is a constant of order 0(Ax/c) when 
y is the fine mesh direction. 


Results 


Representative calculations of invis id and 
viscous flows using the second-order mecht , 

UW2II, are shown in Figs. 1-9. These results are 
compared with chose obtained from other methods 
and exp er imen c s . 

Results for che inviscid transonic flow about 
an 18Z-thick, circular-arc airfoil are illustrated 
in Fig. 1, where surface-pressure distributions 
using the present method and method B2 of Ref. 1 
are given. The two methods show good agreement 
except in che vicinity of che shock wave, where che 
present method shows a superior, three-point 
capture. 

Computed surface pressures for the inviscid 
transonic flow about an NACA 0012 airfoil are 9hown 
tn Fig. 2. This figure shows results obtained usin'* 
che present method and chose obtained from Ref. 2, 
using the central differencing and the (second- 
order upwind) flux-splitting methods. The three 
methods are in basic agreement except in Che neigh- 
borhood of the shock wave, where the present method 
shows a higher quality, two-point capture. 


Utilizing upwind differencing for the inviscid 
term: A3 X , and^central differencing for the 

viscous term v3 , che diagonal difference opera- 
tor tn this expression can be represented by 

(I + (At/2)(A(V x + 4 x ) - ( | A | Ax + 2v)A 7 ]} (25) 


Figure 3 shows calculations of a shock- 
reflection problem investigated by Yee et al. 7 
Figures 3a and 3b show contour plots and pressure 
distributions using the present method UW2II. 
Figures 3c and 3d show results using che present 
method with che high-resolution option added, 
UW2IIHR. (This is the only case in which che high- 
resolution option was used) . These results may be 
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compared with the results of Ref. 12, which is con- 
tained in the present volume, and those of Ref. 7. 
It may be stated that the present method, with the 
high-resolution option off, gives essentially the 
same results as those obtained using the second- 
order flux-splitting method. 7 This is related to 
the fact that the two methods reduce to a single 
method in the linear constant-coefficient case. 

The high- resolution option of the present method 
gives results roughly comparable to those of 
Hartens method (explicit TVD) and to those of the 
implicit TVD method (FTVD) of Yee et al. 12 An 
advantage of the present implicit method and the 
implicit TVD method over the explicit TVD method 
is that the implicit and implicit TVD methods do 
not suffer an explicit time-step limitation and can 
be run at Courant numbers greater than 1. The 
present method was run at a Courant number of 5, 
and convergence was achieved in approximately 
60 time-steps. This is comparable to the perfor- 
mance of the implicit TVD method. 

Figure 4 compares calculations of a Laminar 
supersonic flow over a flat plate at a low Reynolds 
number. In this case, a strong interaction at the 
place leading edge results in a shock wave and a 
nonuniform pressure distribution. The present 
method is compared with results from MacCormack’s 
explicit- implicit method 11 (those results were 
provided by J. Viegas, Ames Research Center). 
Although not shown, comparisons of surface-pressure 
and skin-friction distributions indicate that the 
two methods are in close agreement. Figure 4b 
shows pressure distributions normal to the plate, 
intersecting Che shock wave. In this case, the 
shock resolution of MacCormack’s method is degraded 
by comparison with the present method. 

Figures 5-9 show results of transonic-airfoil 
calculations, at high Reynolds numbers, requiring 
the use of turbulence models. Figure 5 shows the 
mesh system used. The mesh is analytically gener- 
ated and consists of a 120 * 50 cell C-grid. A 
sheared conformal mapping is used over the front 
half of the airfoil, and an H-cype mesh is used 
over the back half and in the wake. The mesh 
points are exponentially spaced away from the air- 
foil, and the cell spacing adjacent to the surface 
corresponds to a y+ of approximately 0.8. About 
23 mesh points are contained in the boundary layer 
at the midpoint, and the outer boundary is placed 
at 10 chord lengths from the airfoil. 

Boundary conditions used in the airfoil calcu- 
lations were as follows: no-slip velocity condi- 

tions at the airfoil surface, constant total 
enthalpy and entropy (or constant total pressure) 
conditions around the outer "C" part of the mesh 
boundary, and constant (free-stream) static pres- 
sure conditions along the vertical back boundary. 
Further details on these boundary conditions are 
given in Ref. 1. 

Figure 6 shows results for an NACA 0012 air- 
foil at the same Mach number and angle of attack 
as used in the inviscid results (see Fig. 2). 

Figure 6a shows computed Mach contours, and Fig. 6b 
shows computed pressure distributions. These are 
compared with preliminary experimental results 
(from the Ames high Reynolds Number Channel II 
facility) obtained by J. McDevitt, L. Hand, and 
A. Okuno. The calculations for this case were 
done with the full Navier-S tokes equations, using 
the zero-equation turbulence model of Cebeci and 


Smith, without pressure-gradient correction. In 
the experiment, the wind-tunnel walls were adjusted 
to conform to free-air streamlines so that inter- 
ference effects should be minimized. The (free- 
air) calculations using the present method are in 
good agreement with experiment. Calculations of 
this flow using various two-equation turbulence 
models have also been made and are reported in a 
companion paper. 1 3 

Figures 7 and 8 show results for the RAE 2822 
airfoil 14 which was used as an experimental test 
case at the 1981 Stanford conference on Complex 
Turbulent Flows. 15 Figure 7 corresponds to the 
subcrltical Case 1 of the experiment, and Fig. 8 
corresponds to the supercritical Case 9. Computed 
Mach contours are shown in Figs. 7a and 8a, and 
computed pressure distributions are compared with 
experimental distributions in Figs. 7b and 8b. The 
calculations were done at the nominal test Mach 
number and the geometric angle of attack, and do 
not account for wind-tunnel interference effects. 
They were made using a new two-equation model which 
is reported in Ref. 13. A more detailed comparison 
of computational and experimental results is also 
Included in Ref. 13. 

Computational efficiency of the present 
method is indicated in Fig. 9, which shows residual 
decay and lift convergence histories for the sub- 
critical RAE airfoil calculation. The calculations 
used a spatially varying time-step in the inviscid 
flow 1 and a constant time-step in the viscous 
boundary- layer flow; 3 min were required for 
500 steps on the Cray-1 computer. Since lift con- 
vergence (and drag convergence, which is not 
shown) occurs at about 200 time-steps, the actual 
computing time required for practical applications 
in this case would be 1.2 min. The computer code 
Is partially vectorized (i.e., vector right-hand 
side and scalar left side) and a speed up by a 
factor of 2 may be expected if the complete code 
is vectorized. An additional speed up may be 
achieved by optimizing the local time-step, and 
work in this direction is currently in progress. 


Concluding Remarks 

We have described a class of implicit upwind- 
differencing methods for solving the compressible 
Euler and Navier-S tokes equations. The methods 
are based on the use of local eigenvalues or wave 
propagation speeds to control inviscid spatial 
differencing and are aimed at achieving more accu- 
rate and stable solutions than can be obtained 
using more conventional methods. Of the class of 
methods investigated, one second-order implicit 
upwind method, UV2II, was found to produce accurate 
solutions of transonic-flow problems without the 
need for special treatments at normal shock waves 
or sonic lines. By combining the present method 
with acceleration techniques based on spatially 
varying time-steps, an effective and efficient 
algorithm is produced for solving compressible 
viscous-flow problems. 
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Fig. 2 Comparison of computed surface pressures 
for inviscid transonic flow about an NACA 0012 
airfoil: M ■ 0.75, a * 2°, 80 * 20 cell C-mesh 

(79 * 31 mesh for Ref. 2). 
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METHOD UW2II 


METHOD UW2II HR 



Fig. 3 Comparison of computed pressure contour levels and distributions for inviscid shock-reflection 
problem: M - 2.9, 9 - 29°, 60 x 20 cell H-mesh. a) Contour plots, present method UW2II; b) Pressure 

distributions, method UW2II; c) Contour plots, method UW2II-HR; d) Pressure distributions, method UW2II-HR. 




Fig. 4 Viscous supersonic leading-edge flow, a) Sketch of flow; b) Pressure distributions normal 
to place surface. 
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Fig. 5 Representative (120 * 50) C-mesh system 
used in viscous transonic airfoil calculations. 




a) Computed Mach contours. 


b) Experimental and computed surface pressures. 


Fig. 6 Transonic flow about NACA 0012 airfoil: M ■ 0.75, Re ■ 10 7 , a * 2°, 120 * 50 cell C-mesh, 

Cebeci-Smith O-equation model. 
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Fig. 7 Transonic flow about RAE 2822 airfoil. Case Is M - 0.676, Re • 5.7 * 10 6 , a - 2.4°, 
120 x 50 cell C-mesh, q-v two-equation model. 




a) Computed Mach contours. b) Experimental and computed surface pressures. 

Fig. 8 Transonic flow about RAE 2822 airfoil. Case 9: M - 0.73, Re - 6.5 * 10 6 , CFL - 6, NE - 500, 
TCR - 0.12, a - 3.19, 120 * 50 cell C-mesh, q-w two-equation model. 
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a) 1*2 residual history. b) Lift history. 

Convergence histories for RAE 2822 airfoil computation. Case 1: M - 0.676, Re - 5.7 x 10 6 , 

4°, CFL - 8, 120 x 50 C-mesh, q-v tvo— equation model. 
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